/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.model.femModel.problemDriver;

import java.util.ListIterator;

import br.ufmg.dees.insane.analysisModel.AnalysisModel;
import br.ufmg.dees.insane.materialMedia.degeneration.Degeneration;
import br.ufmg.dees.insane.materialMedia.degeneration.Representation;
import br.ufmg.dees.insane.model.femModel.element.Element;
import br.ufmg.dees.insane.model.femModel.element.ParametricElement;
import br.ufmg.dees.insane.shape.Shape;
import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.INaturalCoords;
import br.ufmg.dees.insane.util.IVector;
import br.ufmg.dees.insane.util.numericalIntegration.GaussPoint;
import br.ufmg.dees.insane.util.numericalIntegration.IntegrationOrder;
import br.ufmg.dees.insane.util.numericalIntegration.IntegrationPoint;

/** 
 * A class representing the technique of selective integration
 * used to prevent the "locking" in the analysis of thin plates
 * of Reissner-Mindlin.
 *
 * @author Samir
 * @since 20/07/2006
 */
public class ReissnerMindlinThinPlate extends ParametricPhysicallyNonLinearSolidMech {

	private static final long serialVersionUID = 1L;

	/* (non-Javadoc)
	 * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getC(br.ufmg.dees.insane.model.femModel.element.Element)
	 */
	public IMatrix getC(Element e) {
	    int numberOfDegreesOfFreedom = e.getNumberOfDegreesOfFreedom();
	    int k = 0;
	    int kCsi = 1;
	    int kEta = 1;
		IMatrix stiffnessMatrixBending = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    stiffnessMatrixBending.setZero();
		IMatrix stiffnessMatrixShear = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    stiffnessMatrixShear.setZero();
		IMatrix stiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    stiffnessMatrix.setZero();
		IMatrix finalStiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    finalStiffnessMatrix.setZero();
		ListIterator gns = ((ParametricElement)e).getDegenerations().listIterator();
	    AnalysisModel analysisModel = e.getAnalysisModel();
	    Shape shape = e.getShape();
	    IMatrix cartesianNodalCoordsMatrix = e.getCartesianNodalCoordsMatrix();
	    IntegrationOrder iO = ((ParametricElement)e).getIntegrationOrder();
	    String element = ((ParametricElement)e).getClass().getName(); 
	    Representation representation = new Representation();
		while(gns.hasNext()) {
			stiffnessMatrix.setZero();
		    Degeneration degeneration = (Degeneration) gns.next();
	        double[] naturalCoordsBending = degeneration.getRepresentation().getNaturalcoords().toDouble();
	        IVector shapeFunction = shape.getShapeFunction(naturalCoordsBending,cartesianNodalCoordsMatrix);
	        IMatrix derivedShapeFunction = shape.getDerivedShapeFunction(naturalCoordsBending,cartesianNodalCoordsMatrix);
	        IMatrix secondDerivedShapeFunction = shape.getSecondDerivedShapeFunction(naturalCoordsBending,cartesianNodalCoordsMatrix);
			IMatrix matrixB = analysisModel.getInternalVariableOperator(derivedShapeFunction, secondDerivedShapeFunction, shapeFunction, cartesianNodalCoordsMatrix);
			IMatrix bBending = this.getInternalVariableOperatorBending(matrixB);
			IMatrix matrixC = degeneration.mountC();
			IMatrix kB= new IMatrix(numberOfDegreesOfFreedom, matrixC.getNumCol());
			kB.mulTransposeLeft(bBending, matrixC);
			stiffnessMatrixBending.mul(kB,bBending);
			double weightBending = degeneration.getRepresentation().getWeight();
			double intFactor = analysisModel.getIntegrationFactor(shapeFunction, cartesianNodalCoordsMatrix, degeneration.getGeometricProperties());
			double jac = analysisModel.getJacobianTransformation(derivedShapeFunction,cartesianNodalCoordsMatrix);
			double coeff = e.getElmCoefficient();
			stiffnessMatrixBending.setScale(weightBending*jac*intFactor*coeff);
			if(k < Math.pow((Math.sqrt(((ParametricElement)e).getDegenerations().size())-1),2)){
			    if(element.equals("br.ufmg.dees.insane.model.femModel.element.Quadrilateral")){
					IntegrationPoint ip = new GaussPoint();
					for(int i = 0; i < kCsi; i++){
						for(int j = 0; j < kEta; j++){
							INaturalCoords nc = new INaturalCoords();
							nc.setCsi1(ip.getCoord((iO.getXiIntegOrder()-2),i));
							nc.setCsi2(ip.getCoord((iO.getXiIntegOrder()-2),j));
							representation.setNaturalcoords(nc);
							representation.setWeight(ip.getWeight((iO.getXiIntegOrder()-2),i)*
							        	 ip.getWeight((iO.getXiIntegOrder()-2),j));
						}
					}
					kEta++;
				}else if (element.equals("br.ufmg.dees.insane.model.femModel.element.Triangular")){
					IntegrationPoint ip = new GaussPoint();
					for(int i = 0; i < kCsi; i++){
						INaturalCoords nc = new INaturalCoords();
						nc.setCsi1(ip.getCoord(iO.getXiIntegOrder()-1,(3*i)));
						nc.setCsi2(ip.getCoord(iO.getXiIntegOrder()-1,(3*i+1)));
						representation.setNaturalcoords(nc);
						representation.setWeight(ip.getWeight(iO.getXiIntegOrder()-1,i));
					}
				}
			    if(kEta > Math.sqrt(((ParametricElement)e).getDegenerations().size())-1){
			    	kEta = 1;
			    	kCsi++;
			    }
				double[] naturalCoordsShear = representation.getNaturalcoords().toDouble();
		        shapeFunction = shape.getShapeFunction(naturalCoordsShear,cartesianNodalCoordsMatrix);
		        derivedShapeFunction = shape.getDerivedShapeFunction(naturalCoordsShear,cartesianNodalCoordsMatrix);
		        secondDerivedShapeFunction = shape.getSecondDerivedShapeFunction(naturalCoordsShear,cartesianNodalCoordsMatrix);
				matrixB = analysisModel.getInternalVariableOperator(derivedShapeFunction, secondDerivedShapeFunction, shapeFunction, cartesianNodalCoordsMatrix);
				IMatrix bShear = this.getInternalVariableOperatorShear(matrixB);
				IMatrix kS= new IMatrix(numberOfDegreesOfFreedom, matrixC.getNumCol());
				kS.mulTransposeLeft(bShear, matrixC);
				stiffnessMatrixShear.mul(kS,bShear);
				double weightShear = representation.getWeight();
				stiffnessMatrixShear.setScale(weightShear*jac*intFactor*coeff);
				stiffnessMatrix.add(stiffnessMatrixShear);
			}
			stiffnessMatrix.add(stiffnessMatrixBending);
			IMatrix matrixT = analysisModel.getTransformationMatrix(derivedShapeFunction,cartesianNodalCoordsMatrix);
			stiffnessMatrix.mul(matrixT);
			matrixT.transpose();
			matrixT.mul(stiffnessMatrix);
			finalStiffnessMatrix.add(matrixT);
			k++;
			}
		return(finalStiffnessMatrix);
	}
    
	/**
	 * Mounts the internal variable operator bending matrix (BB).
	 * @param B matrix complete 5 x 5.
	 * @return bBending matrix, part of B matrix complete relative the bending.
	 */
    private IMatrix getInternalVariableOperatorBending(IMatrix B){
    	IMatrix bBending = new IMatrix(B.getNumRow(),B.getNumCol());
    	bBending.setZero();
    	
		for(int i = 0; i < B.getNumCol(); i++){
			bBending.setElement(0,i,B.getElement(0,i));
			bBending.setElement(1,i,B.getElement(1,i));
			bBending.setElement(2,i,B.getElement(2,i));
		}
		return bBending;
    }
    
	/**
	 * Mounts the internal variable operator shear matrix (BS).
	 * @param B matrix complete 5 x 5.
	 * @return bShear matrix, part of B matrix complete relative the shear.
	 */
    private IMatrix getInternalVariableOperatorShear(IMatrix B){
    	IMatrix bShear = new IMatrix(B.getNumRow(),B.getNumCol());
    	bShear.setZero();
    	
		for(int i = 0; i < B.getNumCol(); i++){
			bShear.setElement(3,i,B.getElement(3,i));
			bShear.setElement(4,i,B.getElement(4,i));
		}
		return bShear;
    }

}